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Abstract. Motivated by a recent experimental observation of a complex magnetic 
structure [Takada et al. 2013 J. Magn. Magn. Mater. 329 95] we present a theoretical 
study of the magnetic structure of an Fe monolayer deposited on Rh(OOl). We use a 
classical spin Hamiltonian with parameters obtained from ab initio calculations and 
go beyond the usual anisotropic Heisenberg model by including isotropic biquadratic 
interactions. Zero-temperature Landau-Lifshitz-Gilbert spin dynamics simulations 
lead to a complex collinear spin configuration that, however, contradicts experimental 
finding. We thus conclude that higher order multi-spin interactions are likely needed 
to account for the magnetic ordering of the system. 


1. Introduction 

Magnetism on the nanoscale has become a central aspect of modern technology, with 
nanostrnctnres being employed in a mnltitnde of indnstrial applications. For the ongoing 
development of technology and onr basic nnderstanding of the behavionr of novel 
materials experimental and theoretical stndies have to be nsed in concert. Advancements 
in experimental techniqnes, in particnlar, the development of spin-polarized scanning 
tnnneling microscopy (SP-STM) [1], has allowed the direct observation of magnetic 
strnctnre of varions nanosystems at the atomic scale. 

Ab initio calcnlations can greatly help nnderstand the nnderlying physics behind 
observed phenomena, and can hint towards hitherto nnnoticed featnres or direct 
attention to specihc systems. By growing Fe thin films epitaxially on varions transition 
metal snrfaces the effective Fe lattice constant and the hybridization with the snbstrate 
can be tnned, providing deep insight into the relationship between electronic and 
magnetic strnctnre [2, 3]. Using Ir or Rh as snbstrate [4, 5, 6, 7, 8, 9] is especially 
interesting, as their in-plane lattice constant is in between that of bcc-Fe and fcc-Fe. 
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Recent experimental findings by Takada et al. [9] employing SP-STM measnrements 
show that the gronnd state spin conhgnration of the Fei/Rh(001) thin him is a complex 
noncollinear structure with a 4 x 3 magnetic unit cell. A closer inspection of the 
spin structure proposed by Takada et al. [9] reveals that it is not a conventional spin 
spiral, since the spins follow a fanning path rather than a spiral along crystallographic 
directions. 

Earlier theoretical investigations of the same system were based on total energy 
calculations, and considered only a few simple magnetic conhgurations [7, 8]. Ab initio 
methods demand huge computational effort and many aspects of noncollinear magnetism 
as well as hnite temperature behaviour are too difficult to manage. By matching hrst 
principles calculations to classical spin models we may give a less fundamental, yet more 
tractable and more hexible description of magnetic systems. 

In this work we investigate the Fei/Rh(001) thin him in terms of a spin model 
consisting of tensorial Heisenberg couplings, as well as higher order two-spin interactions, 
namely isotropic biquadratic terms. This research is a natural follow-up to our earlier 
study of the closely related Fe/Ir(001) thin hlms [6]. We use the spin-cluster expansion 
(SCE) combined with the relativistic disordered local moment (RDLM) theory [10] to 
obtain model parameters. We also study the ehect of layer relaxations on the magnetic 
interactions. The theoretical spin conhgnration is derived by solving the Landau- 
Lifshitz-Gilbert equations at zero temperature. The highly collinear spin structure 
obtained due to very strong biquadratic interactions allows us to conclude that higher 
order multi-spin interactions are needed to describe properly the magnetism of this 
system. 


2. Methods and computational details 


Computations to obtain the magnetic structure of itinerant systems from hrst principles 
usually involve the adiabatic approximation, assuming a separation of time scales 
between fast (electronic) and slow (spin) degrees of freedom. In terms of the rigid 
spin approximation the orientational state of the spin system is specihed by a set {e} 
of unit vectors e* describing the orientation of the magnetization at site i. 

We study magnetic thin him systems in terms of an extended tensorial Heisenberg 
model of the form 

^({e}) = E 

i i^j 


The traceless symmetric matrices describe the on-site anisotropy energy, whereas 
the exchange tensors F. can be decomposed as 


J.. = 4/ + + J^. , 

—IJ T'J— —IJ —IJ ’ 


( 2 ) 


where / denotes the unit matrix. These contributions can be easily interpreted, since 4 
is the isotropic Heisenberg interaction (note that according to Equation (1) a positive 
Heisenberg coupling is ferromagnetic, i.e., it favours parallel alignment of the interacting 
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spins), the antisymmetric matrix J_^. = ^ {^F. — corresponds to the Dzyaloshinskii- 

Moriya (DM) interaction [12, 13], and J®. = | {^F. + | (^Tr ■ l}j describes the 

second order two-site anisotropy. The last term in Eqnation (1) describes the isotropic 
biqnadratic interaction, with a positive conpling preferring collinear (either parallel or 
antiparallel) orientation. 

We nse the spin-cluster expansion (SCE) introduced originally by Drautz and 
Fahnle [14, 15] combined with the relativistic disordered local moment (RDLM) method 
[16, 17, 18] implemented within the framework of the screened Korringa-Kohn-Rostoker 
(SKKR) multiple scattering theory [11] to obtain parameters for the spin model from 
hrst principles. The SCE allows one to systematically expand the adiabatic magnetic 
energy surface of a spin system over a set of basis functions constructed from spherical 
harmonics, leading to a generalized spin model of the form [14, 15] 

W({e}) = W„ + i: y J^Ydei) 

i L/(0,0) 

+ 5 E E E ly'n(e.)Fi.(ej) + ..., (3) 

^ L^(0,0) iV(0,0) 


where W are real spherical harmonics and the summations exclude the constant spherical 
harmonic of composite quantum number L = {£,m) = (0,0). 

On the one hand, the coefficients of the SCE spin model in Equation (3) can be 
related to those of the conventional spin Hamiltonian in Equation (1). For instance, the 


= 1 components of can be directly related to the tensor, and the isotropic 


biquadratic coupling can be expressed as [6] 

3 




Bii = —- 


Stt 


E J. 

m =—2 


ij 


( 4 ) 


On the other hand, orthonormality of the spherical harmonics implies that the function 
dehned by Equation (3) can be projected to obtain the SCE coefficients, e.g. 


= / d% / d^e,{n)Y,{e,)YL,ie,), 


( 5 ) 


where ('H)g, ^ denotes the partial average of "Hde}) with hxed spin directions at site 
i and j. Associating these averages of the spin model with the partially averaged ab 
initio grand potential allows us to extract interaction parameters from the electronic 
structure. The SCE method is especially useful combined with RDLM theory, wherein 
partial averages of the grand potential can be directly calculated as opposed to using a 
huge number of individual electronic structure calculations to obtain the averages [6, 10]. 
While the method can be formulated to include multi-spin interactions with ease, the 
computation of such terms is currently beyond the capabilities of our computer code. 

Once a set of model parameters is obtained, these may be used to determine the 
ground state magnetic conhguration. To assess the magnetic ordering we use zero- 
temperature Landau-Lifshitz-Gilbert spin dynamics simulations as a means for energy 
minimization. The various terms of the spin model can be easily manipulated in these 
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simulations, in particular to examine how the biquadratic couplings affect the ground 
state spin structure. 

In our calculations the geometry consisted of a single layer of Fe on a semi-infinite 
Rh fcc-(OOl) substrate with a semi-infinite vacuum region on top, with an in-plane 
lattice constant of 2.6898 A corresponding to the experimental value for bulk Rh. Layer 
relaxations were taken into account by varying the distance of the Fe layer from the 
substrate. The local spin density approximation (LSDA) of density functional theory 
was used according to the parametrization of Vosko, Wilk and Nusair [19]. The atomic 
sphere approximation was employed with an angular momentum cutoff of t'max = 2. For 
the computation of interaction parameters Brillouin zone (BZ) integrations were carried 
out with up to 2485 points in the irreducible wedge of the BZ to ensure numerical 
precision. In spin dynamics simulations the initial state was random, and lattices 
consisting of 64 x 64 sites were used with free boundary conditions, or 128 x 128 in 
cases where the wavelength of the expected ground state (based on the Heisenberg 
interaction) was comparable to the smaller system size. 

3. Results 

3.1. Spin model parameters 

We performed self-consistent-field calculations for the Fe monolayer on Rh(OOl) for 
several values of the (inward) layer relaxation of the Fe layer between 0% and —15%, 
—9% being the experimental value [20]. The Fe local spin moment changes moderately 
from 2.93 /xb to 2.84 /ib in this relaxation range, with 2.87 /ib for the experimental 
geometry. We found that in every case the on-site anisotropy is much smaller than the 
two-site contribution, and that the Fe-Fe interaction parameters show strong dependence 
on the layer relaxation, similar to our previous findings in Fei/Ir(001) [6]. Moreover, 
due to the smaller atomic number of Rh with respect to Ir, the relativistic interaction 
terms in the present system are smaller relative to the isotropic Heisenberg couplings 
than in the case of Ir substrate. 

For a quantitative comparison of energy scales, the dominant value of each 
component of the Heisenberg tensor is collected in Table 1 for a few Fe layer relaxations. 
For 0%, —5% and —9% relaxation the components of the first nearest neighbour (NN), 
while for —15% relaxation those of the second NN exchange tensor are shown, with the 
exception of the DM interaction magnitudes which are maximal for second NN’s in every 
case. The isotropic Heisenberg interaction clearly dominates among the bilinear terms, 
while by increasing the inward relaxation the DM interaction significantly increases. It 
should be noted, however, that while at smaller layer relaxations the isotropic terms are 
much larger than the DM interaction, for the experimental geometry and beyond there 
is only a factor of around 7 between the two contributions, implying that relativistic 
interactions have to be taken into account for a proper description of the system. 

In Figure 1 the isotropic Heisenberg interaction, Jh, and the isotropic biquadratic 
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Table 1. Comparison of the dominant energy scales of the bilinear interactions in 
Fei/Rh(001) for a few values of the Fe layer relaxation, obtained with the SCE method 
in a DLM reference state. J\ D and Jzz — Jxx stand for the isotropic Heisenberg 
interaction, the magnitude of the DM vector and the two-site anisotropy, respectively. 
Note the difference in units used. 


relaxation 

Ji [mRy] 

D [^Ry] 

Jzz - Jxx [itRy] 

0% 

1.90 

36.8 

8.73 

-5% 

1.17 

63.8 

8.04 

-9% 

0.591 

77.6 

7.76 

-15% 

-0.555 

79.4 

6.44 


interaction, Bij^ are plotted versus the interatomic distance for various values of the 
Fe layer relaxation. For ideal (i.e. unrelaxed) geometry there is a strong ferromagnetic 
(FM) Heisenberg coupling between first NN’s, with weakly antiferromagnetic (AFM) 
second and third NN couplings. However, as the Fe inward relaxation is increased, an 
AFM tendency arises in the Heisenberg interactions. The first NN couplings decrease 
in value and become AFM at high relaxations, while the second and third NN AFM 
couplings become gradually stronger. In particular at the experimental geometry (—9% 
Fe relaxation) we can see a weakened first NN FM coupling competing with second and 
third NN AFM couplings of the same magnitude, leading to strong frustration of the 
Heisenberg terms. The same qualitative behaviour was also found in Fei/Ir(001) [6], 
but the AFM tendency is weaker in the present system. 




Figure 1. Calculated isotropic Heisenberg (jh) and biquadratic (By) couplings in 
geometrically unrelaxed and relaxed Fei/Rh(001) thin films versus interatomic distance 
(d) in units of the in-plane lattice constant (a 2 d)- 

The biquadratic interactions are very strongly localized, only the hrst nearest 
neighbour couplings are significant. The first NN values are furthermore only weakly 
sensitive to the layer relaxation, and their positive value implies that they prefer a 
collinear arrangement of interacting spins. For small relaxations the Heisenberg terms 
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most likely prefer FM ordering of the spins, which also complies with the weaker 
biquadratic terms. However, due to the comparative insensitivity of the biquadratic 
terms to relaxation, the bilinear and biquadratic couplings are of similar strength near 
the experimental relaxation. Considering that the spin spiral states preferred by the 
frustrated Heisenberg terms are incompatible with the biquadratic couplings of positive 
sign, we anticipate a very strong influence of biquadratic terms on the ground state spin 
structure due to this additional frustration, to be verified by spin dynamics simulations. 

3.2. Spin dynamics simulations 

For ideal geometry the ground state turned out to be indeed FM, with no frustration 
between the exchange interactions and the biquadratic terms. As the Fe layer relaxation 
is increased, however, the competition between these two types of interactions and the 
frustration of the Heisenberg couplings becomes pronounced, and the ground state is 
no longer FM. The change from FM ground state appears quite suddenly. At —7% 
relaxation the ground state is still FM, while at —8% the spin structure becomes 
a seemingly random collinear ensemble of up and down spins, indicating that the 
biquadratic coupling dominates the interaction landscape. Figure 1 suggests that since 
the magnitude of the biquadratic terms is weakly sensitive to Fe layer relaxation, this 
rapid onset of a new type of ground state is due to the decrease of the first NN FM and 
the increase of second and third NN AFM Heisenberg couplings, leading to frustration 
of the exchange interactions. 

This qualitative picture holds true for larger relaxations including the experimental 
geometry. The exchange interactions by themselves would prefer some kind of spin 
spiral due to competing FM and AFM couplings (see Figure 1). This frustration allows 
the comparatively strong first NN biquadratic couplings to overcome the bilinear terms, 
leading to the collinear spin structure shown in Figure 2 for experimental geometry, 
as the biquadratic couplings with positive sign prefer a collinear, either parallel or 
antiparallel, orientation of interacting spins. 

The biquadratic terms themselves do not distinguish between parallel and 
antiparallel pairs of spins. Omitting bilinear terms, the collinear ground state would 
be massively degenerate even for a given common axis of orientation, with spins 
pointing randomly along the common axis. Even though the Heisenberg terms cannot 
destroy the collinear spin structure preferred by the biquadratic terms, they can lift the 
aforementioned degeneracy, appearing as a periodic modulation in the random collinear 
state of spins. 

To trace this effect we computed the lattice Fourier transform of the spin structures 
obtained from simulation, 

1 ^ 

( 6 ) 

' i=i 

where N is the number of spins in the sample and Sj is the unit vector describing the 
direction of the spin at site Rj. The scalar quantity m(q) = ■ m{q) then plays 
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Figure 2. Approximate ground state spin configuration simulated with spin-model 
parameters calculated for experimental geometry. Spins are coloured according to z 
component, with red colour pointing upward and blue colour pointing downward. 


the role of an indicator of any modnlation with wave vector q in the simnlation. In 
Figure 3(a) this scalar is plotted in the entire BZ for experimental geometry. While the 
value of m{q) is overall very small in the sample indicating that there are only weak 
correlations in the structure, there is a pattern emerging from the near-zero background 
in the shape of a rotated square. This is the additional modulation arising from the 
Heisenberg interaction buried in the spin conhguration. 



Figure 3. (a) The scalar lattice Fourier transform, m{q), of the obtained spin structure 
for experimental geometry, (b) The maximal eigenvalues, J{q), of the lattice Fourier 
transform of the calculated exchange tensors for experimental geometry in mRy units. 


To link the pattern seen in Figure 3(a) to the bilinear couplings we determined the 
spatial modulation preferred by them. To this end, we used a mean held estimate based 
on the x(g) paramagnetic spin susceptibility, given by 




SksTl- g{q) 


(7) 
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where T is the temperature, ks is the Boltzmann constant and J_{q) is the lattice Fourier 
transform of the exchange tensor [6]. This formula implies that the highest temperature 
where the susceptibility is singular and, thus, a mean field estimate for the wave vector 
of the ordered state to which the paramagnetic state is unstable, is given by the global 
maximum of the eigenvalues of J_{q) in the BZ. The maximal eigenvalues, J{q), of the 
J_{q) matrices for experimental geometry are plotted in Figure 3(b). 

The similarity between the two quantities shown in Figure 3 is evident. The J{q) 
shows a shallow, nearly degenerate maximum line along a rounded circle, indicating 
that the ground state preferred by the Heisenberg terms is somewhere along these q 
points. It is clear that this maximum line fits perfectly to the Fourier transform of the 
spin stucture in Figure 3(a), verifying that the random collinear configuration forced by 
the biquadratic terms is further modulated by the bilinear interactions. The qualitative 
picture seen in Figure 3 is the same for every relaxation equal to or larger than —8%. 
The sensitivity of the Heisenberg terms to layer relaxation causes the J{q) surfaces to 
evolve somewhat with increasing inward relaxation, leading to a change in the weak 
modulation of the m{q) pattern. 

We also performed a set of spin dynamics simulations with the biquadratic terms 
artificially turned off. As expected, for the larger layer relaxation values the simulations 
evolved into clear single-g helical spin spirals propagating along the (1,1) direction, 
with wave vectors agreeing neatly with the numerical maxima of the corresponding J(q) 
surfaces, being q = (0.47, 0.47) ^ for experimental geometry. The frustration arising 
from competing first and second NN Heisenberg interactions is overall quite similar to 
what we found for the Fei/Ir(001) system [6]. 

4. Discussion and conclusions 

In summary, we found that the (bilinear) exchange interactions in Fei/Rh(001) depend 
strongly on the distance between the Fe overlayer and the substrate, and that their 
competition leads to strong frustration near the experimental geometry, similarly to 
what was found in Fei/Ir(001) [5, 6]. More importantly, we found that biquadratic 
couplings are comparable to the Heisenberg terms, and even dominate those near 
the experimental geometry. Consequently, for Fe layer relaxations between —8% and 
—15% we obtained a complex collinear ground state spin structure from spin-dynamics 
simulations. The fact that the ‘correction’ of the biquadratic terms to the second-order 
spin model drastically alters the ground state spin configuration implies that the usual 
anisotropic Heisenberg model is insufficient to describe this system, and one should not 
even expect to reproduce experimental findings using this approximation. 

It is known that multi-spin interactions in small magnetic clusters can be at an 
energy scale comparable with the bilinear couplings [21]. Their inclusion into the 
mapping of ab initio total energy might significantly affect the computational estimates 
for Heisenberg interactions [22]. If strong enough, these terms can even influence the 
ground state spin configuration of magnetic monolayers [2, 23], in extreme cases leading 
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to the formation of exotic noncollinear spin structures such as nanoskyrmion lattices 

[24]. 

From the point of view of the SCE, a complete second-order (four-spin) SU(2)- 
invariant correction to the Heisenberg model is given by a general energy term [15] 

^(2) = E + E 4jk (ei ■ ej) {ej ■ e^) 

ij ij k 

+ E JISi (e* ■ ej) (Cfc ■ ei ), (8) 

ijkl 

where the summations always run over sets of distinct sites. This would suggest that the 
biquadratic terms we chose to be included in Equation (1) are of the same importance 
as the three- and four-spin interactions. However, in the framework of the SCE-RDLM 
method two-site interactions describe two electron propagations to lowest order, while 
three- and four-spin interactions need at least three and four propagations, respectively 
(cf. Equations (26) and (27) in Reference [10]). Due to the decay of the Green’s 
function with distance, this suggests that the biquadratic terms should, in general, 
be more important than multi-spin interactions, albeit the latter can easily be of the 
same magnitude as multiple-scattering corrections to the former. 

All things considered, the spin structures obtained with our simulations are still very 
unlikely especially with regard to experiments, suggesting that even with the inclusion of 
the isotropic biquadratic terms we are far from grasping the true magnetic ground state 
of the system. As for the spin conhguration observed in the experiment, upon closer 
examination one may realize that the 4x3 period in Reference [9] does not correspond 
to a single-g spin spiral. A lattice Fourier transform of the proposed conhguration seems 
to show modulations corresponding to three distinct, symmetry unrelated wave vectors, 
namely q = (1,0) (^0, ^ and (^4^ 4^ This fact is in accordance with our 

Endings, in that it seems improbable to be dehned by the ground state of a simple 
Heisenberg model. 

5. Acknowledgments 

AD wishes to thank Leonid Sandratskii for an eye-opening discussion. Financial support 
was provided in part by the European Union under FP7 Contract No. NMP3-SL-2012- 
281043 FEMTOSPIN (AD, KP, LS). The work was also supported by the projects 
TAMOP-4.2.4.A/1-11-1-2012-0001 (AD) and TAMOP-4.2.2.A-11/1/KONV-2012-0036 
(LS, IAS) co-hnanced by the European Union and the European Social Fund. 

References 

[1] Wiesendanger R 2009 Rev. Mod. Phys. 81 1495 

[2] Hardrat B, Al-Zubi A, Ferriani P, Bliigel S, Bihlmayer G and Heinze S 2009 Phys. Rev. B 79 

094411 

[3] Simon E, Palotas K, Ujfalussy B, Deak A, Stocks G M and Szunyogh L 2014 J. Phys.: Condens. 

Matter 26 186001 


Magnetic correlations beyond the Heisenberg model in an Fe monolayer on Rh(001) 10 


[4] Martin V, Meyer W, Giovanardi C, Hammer L, Heinz K, Tian Z, Sander D and Kirschner J 2007 

Phys. Rev. B 76 205418 

[5] Kudrnovsky J, Maca F, Turek I and Redinger J 2009 Phys. Rev. B 80 064405 

[6] Deak A, Szunyogh L and Ujfalussy B 2011 Phys. Rev. B 84 224413 

[7] Spisak D, Hafner J 2006 Phys. Rev. B 73 155428 

[8] Al-Zubi A, Bihlmayer G and Bliigel S 2011 Phys. Rev. B 83 024407 

[9] Takada M, Gastelois P L, Przybylski M and Kirschner J 2013 J. Magn. Magn. Mater. 329 95 

[10] Szunyogh L, Udvardi L, Jackson J, Nowak U and Chantrell R 2011 Phys. Rev. B 83 024401 

[11] Zabloudil J, Hammerling R, Szunyogh L and Weinberger P 2005 Electron Scattering in Solid Matter 

(Springer Series in Solid State Sciences, Vol. 147) (Heidelberg: Springer). 

[12] Dzyaloshinskii I 1958 J. Phys. Chem. Solids 4 241 

[13] Moriya T 1960 Phys. Rev. 120 91 

[14] Drautz R and Fahnle M 2004 Phys. Rev. B 69 104404 

[15] Drautz R and Fahnle M 2005 Phys. Rev. B 72 212405 

[16] Gyorffy B L, Pindor A J, Staunton J B, Stocks G M and Winter H 1985 J. Phys. F: Met. Phys. 

15 1337 

[17] Staunton J B, Ostanin S, Razee S S A, Gyorffy B L, Szunyogh L, Ginatempo B and Bruno E 2004 

Phys. Rev. Lett. 93 257204 

[18] Staunton J B, Szunyogh L, Buruzs A, Gyorffy B L, Ostanin S and Udvardi L 2006 Phys. Rev. B 

74 144411 

[19] Vosko S H, Wilk L and Nusair M 1980 Can. J. Phys. 58 1200 

[20] Begley A M, Kim S K, Jona F and Marcus P M 1993 Phys. Rev. B 48 1786 

[21] Antal A, Udvardi L, Ujfalussy B, Lazarovits B, Szunyogh L and Weinberger P 2007 J. Magn. 

Magn. Mater. 316 118 

[22] Lounis S and Dederichs P H 2010 Phys. Rev. B 82 180404(R) 

[23] Kurz Ph, Bihlmayer G, Hirai K and Bliigel S 2001 Phys. Rev. Lett. 86 1106 

[24] Heinze S, von Bergmann K, Menzel M, Brede J, Kubetzka A, Wiesendanger R, Bihlmayer G and 

Bliigel S 2011 Nature Physics 7 713 



